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Topography modeling in atmospheric flows using 
the immersed boundary method 

By I. Senocak, A.S. Ackerman j, D.E. Stevens J and N.N. Mansour 


1. Motivation and objectives 

Numerical simulation of flow over complex geometry needs accurate and efficient com- 
putational methods. Different techniques are available to handle complex geometry. The 
unstructured grid and the multi-block body-fitted grid techniques have been widely 
adopted for complex geometry in engineering applications. In atmospheric applications, 
terrain fitted single grid techniques have found common use. Although these are very ef- 
fective techniques, their implementation, coupling with the flow algorithm, and efficient 
parallelization of the complete method are more involved than a Cartesian grid method. 
Oftentimes, the grid generation can be tedious and one needs to pay special attention in 
numerics to handle skewed cells for conservation properties. Researchers have long sought 
for alternative methods to ease the effort involved in simulating flow over complex ge- 
ometry. A good example is the work by Peskin (1977). He has developed the immersed 
boundary method (IBM) to simulate blood flow in a heart/mitral valve system, where 
the boundary is represented by a body force and the equations are solved on a Cartesian 
grid. However, formulating a suitable body force term has proved to be a challenging 
issue. 

In recent years, IBM has been significantly improved. Mohd-Yusof (1997) has pro- 
posed the direct forcing method in which the body force is implicitly taken into account 
by reconstructing the velocity field around the immersed boundary. Essentially, this new 
approach has eliminated the issue of explicit formulation of a suitable body force repre- 
senting the boundary. IBM with the direct forcing approach has been further developed 
by adopting higher order reconstruction schemes, and it has been successfully applied 
to complex flow problems in engineering applications (e.g. flow over a truck, flow in 
piston cylinder assembly and flow in a stirred tank) by Fadlun et al. (2000); Verzicco 
et al. (2000); Iaccarino & Verzicco (2003). For high Reynolds number flow simulations, 
IBM coupled with an adaptive mesh refinement technique has been effective to provide 
the required near- wall resolution (Kalitzin & Iaccarino 2003) . Adaptive mesh refinement 
helps reduce the overall number of computational nodes to achieve the desired near-wall 
resolution, but the problem can still be very demanding for higher Reynolds number flow 
simulations. 

Wall/surface modeling in LES or the wall-function formulation in Reynolds-averaged 
Navier-Stokes (RANS)computations have been proposed to decrease the computational 
load due to near-wall resolution. In these approaches, coarse resolution grids can be 
employed, and no-slip conditions are not applied directly at the surface, because the 
implied stress would be overestimated on a coarse resolution grid. Instead, stresses at 
the surface are imposed as boundary conditions, alleviating the need to resolve the thin 
turbulent boundary layer. If one wants to extent the IB method for topography modeling 
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in LES of atmospheric flows, for which numerical resolutions have been poor, then it is 
necessary to develop a reconstruction scheme that can take into account surface modeling. 
This issue has not been addressed within the context of IBM and it is one of our main 
objectives in this paper. 

In what follows, we briefly describe the immersed boundary method, and apply it to 
low Reynolds number laminar flow cases to test our implementation. To take into the 
LES surface modeling issue within the IB method, we develop a reconstruction scheme 
based on the mean logarithmic wind profile assumption near the surface. To test this new 
scheme, we perform LES of a neutrally stratified atmosphere. Specifically, we compare 
the results of the IB method with the new reconstruction scheme to the results of the 
commonly adopted surface modeling approach, in which surface stresses are imposed as 
boundary conditions. 


2. Governing equations 

The governing equations for LES of a neutrally stratified atmospheric boundary layer 
are the filtered Navier-Stokes equations written as follows 


duj _ dvn d (ujUj) _ dp_ _ dnj 

dxi ’ dt dxj dxi tljk jUk dxj 


(2.1) 


where turbulent stresses are defined as = uju] — UjUj, fj is the Coriolis parameter, 
tijk is the permutation tensor, u and p are the filtered velocity and filtered dynamic 
pressure, respectively. 

To solve the above governing equations, we use the LES atmospheric research code 
DHARMA. The numerical method adopted in DHARMA is described in detail in Stevens 
& Bretherton (1996) and Stevens et al. (2000). The governing equations are integrated us- 
ing a forward-in-time projection method based on an explicit second-order Runge-Kutta 
scheme (Bell & Marcus 1992). The spatial discretisation is performed on a staggered 
grid. A third-order accurate upwind-biased monotonic scheme is used for the advection 
terms, whereas diffusion and pressure gradient terms are discretized using second-order 
accurate central differencing schemes. A direct solver (FFT) is utilized for solving the 
pressure Poisson equation, and the code has been parallelized to run on various platforms 
using MPI. 


3. Sub grid scale turbulence modeling 

Since the Reynolds number of a typical atmospheric boundary layer (ABL) flow is very 
high, LES of ABL with near-surface resolution is not a practical option with current 
computer resources. To alleviate this obstacle one resorts to LES with surface modeling. 
For a planar surface, a common approach is to set the vertical component of the velocity 
to zero, and to define the horizontal components of the turbulent stresses based on the 
mean logarithmic wind profile assumption (Moeng 1984). The surface stresses can be 
written as follows 

Ti3 = ~ M* <> (3- 1 ) 

where z i, Ui are the distance in the vertical direction and the velocity of the first grid 
point away from the surface, respectively, zo represents the roughness height. It should 
be noted that, in the immersed boundary method, we do not utilize the above approach 
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Figure 1 . Schematics of the reconstruction schemes. Bilinear inverse interpolation (left), 

log-law (right) 


for surface modeling. Instead, we employ the log-law reconstruction scheme, which is 
described in the next section. 

Turbulent stresses that appear in equation 2.1 are modeled based on the Boussinesq 
eddy viscosity assumption. As discussed in Senocak et al. (2004), in LES of ABL near- 
surface models are needed to improve the predictions in the vicinity of the surface. Based 
on the results of that study we employ the hybrid RANS/LES model for the present 
computations. The hybrid RANS/LES model adopts the Prandtl’s mixing length model 
(Prandtl 1925) near the surface and blends in with the dynamic Smagorinsky (Germano 
et al. 1991) away from the surface. The resulting form of the turbulent eddy viscosity 
can be written as 

v t = [(1 - exp (-z/h)) 2 • (CA 2 ) + exp(— z / h) 2 (kz) 2 ]\S\, (3.2) 

where z is the distance from the surface, A is the filter width, k is the von Karman 
constant with a value of 0.41, and |S| is the magnitude of the strain rate tensor. The 
dimensionless parameter C is computed dynamically during the solution, making it a 
function of space and time (Germano et al. 1991). In the above equation h is the altitude, 
corresponding roughly to the upper edge of the surface layer. A value of 150(m) is used 
for the present computations. It should be mentioned that a logarithmic velocity profile 
is expected within the surface layer, which is approximately the bottom 10 percent of 
the atmospheric boundary layer height for neutrally stratified conditions (100 m -200 m) 
(Stull 1988). 


4. Immersed boundary method 

In IBM the solid boundary is represented by a body force F). The discretized form of 
the momentum equation given in 2.1 can be written as follows. 

8 At 1 = RHSi+F i' (4- 1 ) 

where RHSi includes the pressure gradient, convective, diffusive, and Coriolis terms. F) 
is the body source term that gives the desired velocity at the boundary. In the direct 
forcing technique (Mohd-Yusof 1997; Fadlun et al. 2000) if the desired velocity at the 
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boundary is w" +1 = Vfl +1 then one can write the explicit form of the body force as 

V n+1 - w” 

Fi = -RHSi + 8 M V (4.2) 

Hence, instead of imposing the body force F t explicitly to obtain the desired velocity at 
the boundary, one can impose the desired velocity and take into account the body force 
implicitly. For complex geometry, the boundary is not coincident with the Cartesian grid 
nodes, so one needs to reconstruct the velocity field using the values from neighboring 
nodes and the desired value at the boundary. 

The steps involved in applying the immersed boundary method can be summarized as 
follows. 

(a) Preprocessing: Determine the Cartesian cells that are cut by the boundary. Tag 
the nodes as dead , fluid and cut. 

( b ) Predictor stage: Solve the discretized momentum equations. 

(c) Set zero velocity field on the dead nodes, and apply the reconstruction scheme on 
the cut nodes. 

(d) Solve the pressure Poisson equation 

(e) Update the velocity and the pressure field, and impose the reconstruction on the 
cut nodes. 

Several reconstruction schemes have been suggested in literature (Fadlun et al. 2000; 
Iaccarino & Verzicco 2003). The simplest reconstruction scheme is the nearest neighbor. 
For instance, if we consider a 2-D geometry, a linear interpolation can be employed in 
either the x (horizontal) or the z (vertical) direction depending on the distance from the 
boundary. To preserve the local maxima or the minima a bilinear inverse interpolation 
can also be adopted. As depicted in the left part of figure 1, a quadrilateral element can 
be constructed from the neighboring nodes and the boundary surrounding the cut node. 
This quadrilateral element and the coordinate of the cut node is then mapped onto a 
square element, where an area weighted average is used to interpolate the value of the 
velocity at the cut node. A quadratic equation needs to be solved for this mapping, which 
can be derived as follows 

4 4 

xp = ^2xiNi, yp = ^yjNj, 
i—1 i= 1 

Ni = |(1 +££i)(l + VVi)f & = ±Mi = ± 1. (4.3) 

Eliminating £ gives a quadratic equation for r]. 

arf + br] + c = 0, (4.4) 

Upon solution of this equation the root that lies within the mapped zone is taken as the 
solution. A stable algorithm for bilinear inverse interpolation is given by Felippa (2004). 
However, the extension of this method to 3-D (trilinear inverse interpolation) is not 
straightforward. An iterative algorithm is needed to solve the set of equations defining 
the mapping of a hexahedral element to a cubic element 

For high Reynolds number flows, both the nearest neighbor and the bilinear inverse 
interpolation schemes would require a very fine resolution near the boundary. The com- 
putational grid should be fine enough to resolve the viscous sublayer, where the velocity 
variation is physically linear. However, this is an overwhelming requirement, which is not 
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practical for atmospheric flows. Instead one should consider LES of atmospheric flows 
with surface modeling. 

We propose the log-law reconstruction scheme for incorporating surface modeling into 
the immersed boundary method. It should be emphasized that this scheme addresses 
the surface modeling issue through reconstruction of the velocity field, whereas in the 
commonly adopted surface modeling approach turbulent stresses are imposed at the 
surface as boundary conditions. In the following, we explain the formulation of the log- 
law reconstruction scheme. 

The log-law for a rough surface can be written as follows (Panofsky & Dutton 1984) 


U_ 
u * 


-ln(-), 
k z o 


(4.5) 


where z is the distance from the surface measured along the surface normal direction, z 0 
is the roughness height, w* is the friction velocity, and U is the magnitude of the velocity. 
Within the log-layer the friction velocity, «, is constant in the surface normal direction. 
Using this property, one can write the following between two points lying on the same 
surface normal direction. 

U 2 Injzg/zo) . 

Ux ln( Zl /z o y 1 J 

The above formula is based on the magnitudes of the velocity. One needs to decompose it 
into velocity components in each direction (u, v, w), in order to impose the direct forcing 
in equation 4.1. For instance, if we consider a neutrally stratified atmospheric boundary 
layer over a planar surface, the velocity direction in the horizontal plane 9 changes with 
altitude due to the rotation of the Earth. Since one has the magnitude of the velocity at 
hand due to equation 4.6, 9 needs to be extrapolated from the fluid nodes to the cut node 
lying along the same vertical direction. The direction of the velocity on the horizontal 
plane is computed as follows 


6 = arctan(v / u) (4.7) 

For simplicity, we suggest a linear extrapolation to compute 0 at the cut node. Once 9 cu t 
is determined, then horizontal components of the velocity is computed as shown below 

u = U cos 6 cu t , v = U cos 6 cu t (4.8) 

The log-law reconstruction scheme in the case of a planar surface is described schemat- 
ically on the right part of figure 1. In the preprocessing stage, first, the nodes are tagged 
as dead, fluid and cut. Next, the magnitude of the velocity at the cut node, numbered as 
1, is updated based on equation 4.6, using the information from the fluid node numbered 
as 2. Then 0 at fluid nodes 2 and 3 are computed based on equation 4.7. Using 0 2 and 63 , 
a linear extrapolation gives 9\ , which is then used to compute the horizontal components 
of the velocity at the cut node, based on equation 4.8. Note that, the vertical component 
of the velocity (w) has zero value at the boundary and a linear interpolation is used to 
impose it on the cut node. 

For a three dimensional complex surface the above scheme needs to be reformulated 
based on normal and tangential directions to the surface. We are currently working on 
extension of our scheme to complex topography. 



ah 


336 


I. Senocak et al. 



Figure 2. Time evolution of the recirculation zone length and the flow structure. : 

present computations, • : experiment (Coutanceau & Bouard 1977a). 



Figure 3. Streamwise component of the velocity field in the wake of the circular cylinder. 

: present computations, • : computations of Nieuwstadt & Keller (1973), A : experiment 

(Coutanceau & Bouard 1977a). 

5. Results 

In this section we, first, present results of low Reynolds number laminar flow compu- 
tations adopting the linear reconstruction scheme in the IB method. Following this, we 
present results of LES of a neutrally stratified atmospheric boundary layer adopting the 
log-law reconstruction scheme in the IB method. 

5.1. Linear reconstruction scheme 

To test our implementation of the immersed boundary method, we perform simulations 
of laminar flow past circular cylinders and compare our results with the available exper- 
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U/u, 



U/u . 


Figure 4. Comparison of the mean wind profiles. Linear plot (left), semi-log plot (right). 
: surface-stress boundary condition, : immersed boundary with log-law reconstruc- 
tion scheme, : log-law 


imental and computational studies in literature. We consider flows at Reynolds number 
of 20 and 40. The Reynolds number is defined as 


Re = 


UD 


(5.1) 


where U is the upcoming freestream velocity, D is the diameter of the circular cylinder, 
and v is the kinematic viscosity of the fluid. The wake behind the circular cylinder 
becomes unsteady above Reynolds number of approximately 40, (Coutanceau & Bouard 
1977a, 6). 

Figure 2 shows the time evolution of the length of the recirculation zone that forms 
in the wake at Reynolds number of 20 and 40. In figure 3, the streamvise component 
of the velocity is compared with both the experimental data of Coutanceau & Bouard 
(1977a)and the computational data of Nieuwstadt & Keller (1973) .The present results 
obtained with the IB method agrees well with both the experimental and the computa- 
tional data. We have also compared the results of linear reconstruction scheme with the 
results of bilinear inverse reconstruction scheme. Since the results are nearly identical, 
we do not include the comparisons in this study. 


5.2. Log-law reconstruction scheme 

To test the performance of IB method in LES of atmospheric flows, we consider a neutrally 
stratified atmospheric boundary layer. We compare the results of the IB method with log- 
law reconstruction scheme to the results obtained by adopting the surface-stress bound- 
ary condition approach, which is described in equation 3.1. For surface-stress boundary 
condtion approach, the computational domain size is 3000 m x 1500 m x 1500 m with 
64 x 32 x 64 grid points uniformly distributed in x, y, and z directions, respectively. For 
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Figure 5. Comparison of streamwise component of the total (left) and resolved (right) stress. 
: surface-stress boundary condition, : immersed boundary with the log-law recon- 
struction 


immersed boundary simulations, we consider a domain height of 1687.5m with 72 grid 
nodes in the vertical direction. We place the immersed boundary at a height of 187.5m 
so that we have a computational domain, which is identical to the computational domain 
adopting the surface-stress boundary condition approach. The flow is driven by a mean 
pressure gradient that would balance a 10 ms -1 geostrophic wind in the x direction. 
Coriolis parameter is set equal to 10 -4 s -1 , and the roughness parameter zq has a value 
of 0.1 m. A dimensionless time unit can be defined based on the Coriolis parameter(l//). 
The simulations were run over a period of 10/ _1 , and statistics are collected during the 
time period of last the 4/ _1 . Ensemble averaged vertical profiles have been obtained by 
collecting data at every 6 x 10 -3 dimensionless units (60secs), and averaging them both 
in time and in horizontal space. 

In figures 4 and 5 the mean wind profile and the streamwise component of the total 
and the resolved stress obtained with the commonly adopted surface-stress boundary 
condition are compared with the results of IB method. Clearly, the IB method with the 
log-law reconstruction is able to reproduce the results of the surface stress boundary 
condition. 

Figure 6 compares the spanwise components of the total stress and the velocity. The 
noticable differences seen in this plot is because of the extrapolation of the velocity direc- 
tion in the IB method with log-law reconstruction, given in 4.7. As discussed in Senocak 
et al. (2004) the adoption of near-surface models results in non-vanishing spanwise com- 
ponent of the velocity. In that respect, the IB method with log-law reconstruction does 
a better job by predicting a lower value of the spanwise velocity. 

Finally, we compare turbulent eddy viscosity profiles because it is an important quan- 
tity in atmospheric modeling. In figure 7, we see that the IB method with log-law recon- 
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Figure 6. Comparisons of spanwise component of the total stress (left) and the velocity (right). 
: surface stress boundary condition, : immersed boundary with the log-law recon- 
struction 



Figure 7. Comparison of turbulent eddy viscosity profiles. : surface stress boundary 

condition, : immersed boundary with the log-law reconstruction 

struction gives results that are almost identical to results of the surface-stress boundary 
condition. 


6. Summary and conclusions 

In this study, we have presented a Cartesian grid immersed boundary method for LES 
of a neutrally stratified atmospheric boundary layer flow. The IB method is appealing 
for atmospheric flow modeling because its implementation is easy and does not alter the 
computational structure of the Cartesian grid code, which can involve extensive modeling 
for cloud microphysics and radiation modeling. 

For validation purposes, first, we have considered laminar flows over a circular cylinder 
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at low Reynolds numbers and demonstrated good agreements with available experimental 
and computational studies in literature. To address the surface modeling issue in LES of 
ABL using the immersed boundary method, we have proposed the log-law reconstruction 
scheme. This scheme enables us to employ the IB method for high Reynolds number flows 
without the need for fine near-surface resolution. We have shown that the IB method with 
log-law reconstruction scheme can produce results that are nearly identical to the results 
obtained with the commonly adopted surface-stress boundary condition approach. 

Our future work will focus on applying the IB method to atmospheric flow problems 
that involve complex topography, cloud coverage and radiative energy transfer. 
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